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, A numerical simulation of a gas-fluidized bed is performed without introduc- 

tion of any empirical parameters. Realistic bubbles and slugs are observed in 



X 



our simulation. It is found that the convective motion of particles is impor- 
tant for the bubbling phase and there is no convection in the slugging phase. 



From the simulation results, non-Gaussian distributions are found in the par- 
ed 



ticle velocities and the relation between the deviation from Gaussian and the 
local density of particles is suggested. It is also shown that the power spectra 



of particle velocities obey power laws. A brief explanation on the relation- 
ship between the simulation results and the Kolmogorov scaling argument is 
discussed. 
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I. INTRODUCTION 



Recently dynamics of granular systems has attracted much attention among physicists 
[ l]^5|] as a typical object of non-equilibrium statistical physics. For example, in vibrating 
beds |6|-[T0|, collisions among particles produce convection and turbulence. However, the 
fluidization of granular particles immersed in a fluid stream, where the hydrodynamic inter- 
actions are relevant, exhibits richer phenomena, illustrative of dynamical phase transitions 

EMI- 

In an experiment on fluidized beds, we prepare a vessel containing granular particles 
and impose a gas flow from the bottom of the vessel. When the flow rate is small enough, 
particles do not move. This state is known as a fixed bed. Above a critical value of flow 
rate, the fixed bed is destabilized, and then is fluidized uniformly. At larger flow rate, 
the uniformly fluidized bed becomes unstable and bubbles appear. Increasing the flow rate 
further, the bubbles become larger and then become slugs which are horizontally spread 
bubbles. For further increase of flow rate the state becomes disordered, and finally reaches 
a dilute state of particles in which is recovered spatial homogeneity. These phenomena are 
similar to boiling of water. The phase transitions in fluidized beds, however, are not thermal 
phase transitions. Thus the mechanism of the phase transitions in fluidized beds must differ 
from that of the boiling of water. 

These dynamical phase transitions have not been observed in the systems of smaller 
particles in flows, such as colloid particles which are much smaller than granular particle 
|14| . The reason such interesting phenomena are observed in granular systems may be 
absence of Brownian motion. Namely, granular systems cannot reach any equilibrium states, 
while colloid particles can reach. Thus, granular particles are one of interesting subjects in 
nonequilibrium statistical physics. 

We do not have any established model which is suitable to describe fluidized beds. A mod- 
ern and successful approach is simulation by the distinct element method (DEM) |15|-[T7|. 



The DEM is also a powerful tool to describe vibrating beds PJTOfl. In this approach, the 
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interactions among particles are replaced by a mechanical model which consists of springs, 
dashpots and sliders. The fluid motion is assumed to obey a phenomenological model where 
the hydrodynamic interactions are replaced by the coarse-grained friction between the par- 
ticles and the fluid. In this way, hundreds of thousands of particles have been successfully 
simulated and a realistic motion of particles is reproduced. However, we stress that the 
DEM is not independent of the experiments because experimental results are used to choose 
parameters. 

Two-fluid models are often used to describe fluidized beds ||l2l , |l3| . These models allow one 
to understand macroscopic pattern formation, using bifurcation analysis and hydrodynamic 
stability analysis. In addition, simulations based on the two-fluid models reproduce realistic 
motions [|T~3| , p~8| , p~9|] , and some authors indicate that solitons play an important role near 
the onset of the instability of uniformly fluidized beds p0|-f2^]. In spite of these successful 
results, the two-fluid model contains some difficulties. For example, it is difficult to choose 



a suitable two- fluid model [2^] . The role of particle motion is not clear because particles are 
described as a fluid. Furthermore, the two-fluid model is supplemented by empirical laws for 
the choice of parameters. Although two-fluid models contain a continuous approximation, 
the simulation of two-fluid models is not easier than the DEM. 

In this paper, we perform a simulation based on the model which does not contain any 
empirical parameters except for the particle radius, the mass densities of the fluid and the 
particles, and the shear viscosity of the fluid. 

For this purpose, we neglect the complexities of granular systems, which are polydis- 
persity, direct interactions among the particles such as the Coulomb interaction and inter- 
molecular forces, and chemical reactions induced by mixing, etc. We only treat systems 
which contain monodisperse spheres with hard-core interactions among the particles and 
with hydrodynamic interactions described by the Stokes approximation. 

In the next section we show how to construct the model and discuss its relevance in 
detail. In Sec. ^TTJ, it is shown that bubbles and slugs are observed in the simulation. In Sec. 



rV| we analyze the data obtained from the simulation. For example, we show the distribution 



functions and the power spectra of particle velocities. In Sec. [V], we briefly explain how 
power laws in power spectra appear. In Sec. [VT] , we conclude and summarize our results. 
We will present the details of our calculation method for hydrodynamic interactions among 
particles in Appendix A, and the treatment of fixed particles in Appendix B. 



II. SIMULATION METHOD 



In this section, we summarize the algorithm of our simulation for granular particles 
immersed in a fluid stream. This section consists of two parts. The first part is devoted to 
the general perspective of the motion of particles in a fluid. In the next part, we discuss 
the validity of our approximations. The hydrodynamic interactions are calculated by the 
Stokesian dynamics method [f26| -|28H which is briefly described in Appendix 

In general, classical particles with mass m in a fluid obey the Langevin equation 



mjV = F f + F g + F i + F b , 



(1) 



where the velocity U, the position x and the forces F; (7 = f,g,i,b) represent vectors 
containing N particle elements as 
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where the superscript (i) represents the index of the particle. In Eq. (|l]) there are four kinds 
of forces: Ff is the drag force from the fluid, F g is the gravitational force exerted on the 
particles, F, is the force due to direct interactions among particles, and F^ is the Brownian 
force coming from the thermal motion of the fluid. Although the size distribution and the 
shape of particles are important factors in technology, we restrict ourselves to the motion of 
monodisperse spherical suspensions. 

We define the following dimensionless quantities, 
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where Pe,Re and St are respectively the Peclet number, the Reynolds number and the 
Stokes number. In Eqs.(|3])-(|D, o> Pp, Pf and p are the particle radius, the mass densities 
of the particle and the fluid, and the shear viscosity of the fluid respectively, ks is the 
Boltzmann's constant. V is the characteristic velocity, where we choose the sedimentation 
rate of one particle 

u d = ^L = 2 ^M (6) 

birpa yp 

with the effective gravitational acceleration g = g{p p — Pf) / p p . These dimensionless numbers 
actually have definite meanings. The Peclet number is the ratio of work done by the drag 
force over the size of a particle (67ipa 2 V) to thermal energy (/c#T). The Reynolds number 
is well known as the ratio of inertia to drag. The Stokes number represents the relative 
importance of the kinetic energy of particles (mV 2 ) to work done by the drag force. It is 
also recognized as the ratio of the time scales, 

St = % (7) 

where T r = m/Qirpa is the relaxation time of particle velocity due to the drag force, and 
T p = a/V is the passing time of the particle scale a with the velocity V. 

The Froude number, which has been widely used on this problem, is given by 

V 2 

Ft = — , (8) 

gL 

with L being the linear size of system. The meaning of the Froude number is the ratio of 
the kinetic energy to the gravitational potential, and Fr is proportional to St if we adopt 
V = Uq. For later convenience, we also introduce the effective Reynolds number of particles 

as 



where v v = \ij p p . The Stokes number is related to the particle size, while the Froude number 
and the effective Reynolds number Re^ p ) are related to the system size. The role of Re^ p ) will 
be discussed in Sec. |V[ Thus Pe,Re and St are fundamental and independent parameters 
to determine the motion of particles. 

Let us consider the explicit form of forces in (|1|) and discuss their relevance. The gravi- 
tational force F g acting in the z direction is simply given by 

F g = -m~gE z , (10) 

where E 2 is the generalization of unit vector e 2 as in (H). The direct interaction among 
particles Fj is assumed to be due to hard-core interactions. Therefore it can be treated as 
exchanging velocities during elastic collisions among particles. We neglect the random force 
Fb, because in this paper we treat the system where the Peclet number is large enough. The 
validity of Pe ^> 1 will be discussed later. 

The most relevant and complicated force is Ff, which is determined by the Navier-Stokes 
equation 

Pf {d t u+{u- V)u} = ^V 2 u - Vp, (11) 
with the incompressible condition 

V-u = 0. (12) 

The incompressibility is valid even if the fluid is air, when the particle motion is much slower 
than the sound velocity of fluid. When i?e « 1, the Navier-Stokes equation is reduced to 
the Stokes equation 

-/A7 2 u + Vp = 0. (13) 
In this case the particle velocity U is connected with the force induced by the particles on 



the fluid F p in the following linear relation |29 



F p = R (x) • (U - u 00 ) , (14) 

where u°° is the fluid velocity in the absence of particles and R is the 3iV x 3iV resistance 
matrix which depends only on the particle configuration x. The details of the construction 
of R are described in Appendix [A[ The force F p induced by the particle is related to Ff by 

F p = -F,. (15) 

In later discussions, we assume Re < 1. As we will show, this approximation is valid for 
the motion of relatively small particles. 

In the case of Re <C 1 and Pe ^> 1 and with hard-core interactions, the Langevin 
equation (JJ) can be reduced to 

Stittf) = -R(x(t)) • (ti(t) - u°°) - E*. (16) 



dt 

Here we scale the velocities by V, the length by the radius a and the resistance matrix by 
the drag factor 67r/za. We denote the dimensionless value of U as U. From ( fT6| ) it is obvious 
that the Stokes number St is an important parameter for our model. The role of St is also 
pointed out for dilute monodisperse suspensions with low Reynolds number fl3"0| , |3~T |. For 



vibrating beds in which the particle inertia is dominant and the hydrodynamic interaction 
is negligible, St should be large while for liquid-fluidized beds St is small. 

In our simulation, we need to integrate Eq. (|1~6|) for a small time interval numerically. We 
divide U(t) into two parts, 

ti(t) =U + Ui(f). (17) 

Uo is determined by 

= -R(x)-(u Q -u°°)-E z . (18) 

It is obvious that Uo is the same solution as that for St <C 1. Although Ui(t) is determined 
by 
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Stiui® = -R(x) • Ui(t), (19) 

we adopt the simplest form for computational efficiency, 

Ui(t) = Ui(0) expf-^, (20) 

which is the solution of ( jl9|) on the assumption of R ~ I, where I is the unit tensor. From 
the initial condition U(0) = U + Ui(0), we thus obtain 

U(f) = U + (U(0) - U ) exp . (21) 

This solution shows us that the particle velocities are damped from the initial value U(0) to 
the terminal velocity Uo- The simplification in Eq.(|20|) is crucial. We expect, however, the 
error from this simplification is small when we choose a small enough time interval for the 
numerical integration, because the resistance matrix is calculated as a function of particle 
configuration at each step. 

Let us estimate the values of the dimensionless parameters. We adopt p p = 2.5 [g cnr 3 ] 
which is a typical value of glass beads. Fluids are assumed to be air (pf = 1.2 x 10~ 3 [g cm" 3 ] 
and p = 1.82 x 10~ 4 [g/cm sec]) and water (pf = 1.0 [g cm -3 ] and p = 1.0 x 10~ 2 [g/cm sec]) 
at room temperature. Substituting k B = 1.38 x 10~ 16 [erg K" 1 ], T = 293[K] = 20[°C] and 
g = 981 [cm sec -2 ], we obtain Table [I] and Fig. [I]. We also show the characteristic times, T r 
and T p , in ([?]). 

The lines in Fig. [I] connect the points where dimensionless parameters are equal to unity 
for air and water. Therefore, there is no meaning in the vertical axis. In the left of the 
line of Pe = 1 in Fig. [I], t ne random force is dominant and the particles can be regarded 
as typical colloidal suspensions. On the other hand, the particles to the right of Pe = 1 
can be regarded as typical granular particles where the hydrodynamic interaction among 
particles dominates the Brownian force. In the left region of Re = 1, the fluid motion can 
be described by the Stokes equation. On the other hand, in the right of Re = 1, the inertia 
term in the Navier-Stokes equation is important. In the left region of St = 1 where T r < T p , 
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we can neglect the inertia effect of particles. Thus, the particles move with their terminal 
velocities. In the right of St — 1 where T r > T p , the effect of particle inertia exceeds the 
drag from the fluid. Therefore, collisions among the particles are important. 

We choose the particle radius as 10/xm because we are interested in the case of Pe ^> 1 
and Re <C 1. In this case the dimensionless parameters have the values presented in Table 
| for air and water. As anticipated in Fig. [I], the difference between air and water appears 
in the value of the Stokes number. We now interested in the phenomena driven by the 
drag force, which is the behavior in the time scale T p . Therefore, for gas-fluidized beds, the 
collisions among particles are not negligible, while for liquid-fluidized beds, the collisions are 
not important. We focus on gas-fluidized beds in this paper. We will discuss the properties 
of liquid-fluidized beds elsewhere. 

We notice that the particle radius a ~ 10~ 3 [cm] belongs to the group C in the classifica- 
tion by Geldart |32| , pT3| in which all cohesive powders are difficult to fluidize. This difficulty 



in experiments arises because direct interparticle forces dominate the hydrodynamic drag 
force. If the interparticle forces, except for the hard-core interactions can be removed, a 
collection of small particles must be fluidized as will be shown. 
Finally we summarize the adopted assumptions. 

1. the diameter of all particles is identical. 

2. the interaction among particles is described by hard-core interactions. 

3. the thermal motion is negligible. 

4. the inertia of the fluid is negligible. 

Assumptions [5] and [| are valid when the particle radius is between 1 [/im] and 10 [/im] in both 
air and water. 
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III. SIMULATION RESULTS 



In this section, we collect three typical results of simulations. In many of our simulations, 
the configuration of particles is restricted to lie in a plane parallel to the direction of gravity 
for efficiency of calculation, although periodic boundary conditions are also applied in the 
direction perpendicular to this plane. Thus, in such cases, particles are influenced by three- 
dimensional hydrodynamic interactions. 

At first we show the result of our simulation using rectangular cells with periodic bound- 
ary conditions in which the ratio of the height to the width is large compared with unity. 
We observe stable slugs which move upward in this simulation (Fig. §). This result is also 
observed in the three-dimensional (not monolayer) simulation with slender cells as well (Fig. 
|3|). From these results, particles in the dilute region relatively fall down, so that in the 
denser region the sedimentation velocity is smaller. This tendency agrees with the standard 
theory [|33| -|35|| and the experiments of sedimentation |1^J36[| . In our simulation for slugs the 
configuration of particles is nearly closed packed and particles are immobile in concentrated 
regions. On the other hand in the dilute region, particles monotonically fall down. We can- 
not observe any convection of particles in this simulation. Thus slugs can be produced by a 
pseudo one-dimensional motion of particles, after the uniform state has become unstable. 

Next we show the result of our simulation using square unit cells with periodic boundary 
conditions. Figure § is a typical snapshot. In contrast to the previous result, we see that 
particles in the dilute region float up, while particles in the dense region fall down. Thus 
particle convection exists around the dilute region, which may be regarded as a bubble. This 



kind of convection inside bubbles has been observed in experiments p7 |. Thus, we infer that 
convection is important to create bubbles. This bubble, however, is not stable and it will 
disappear soon afterwards. We can also see periodic birth and death processes of bubbles. 

We have also performed simulations introducing fixed particles in our system. The 
treatment of fixed particles is described in Appendix [B]. In real systems, particles are settled 
in a vessel. For this purpose, we introduce particles fixed at the bottom of the unit cell. In 
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our simulation, the fixed particles are placed horizontally, spaced as far as 3.5 radius apart. 
We perform this simulation as follows. At first we randomly position free particles and allow 
them to drift to the bottom under the influence of gravity. Thus the free particles fall down 
and produce a fixed bed. Next we inject the flow upward with the velocity u°°. The result 
of our simulation after the injection of flow is shown in Fig. We notice that the unit cell 
drawn in Fig. || is connected with its mirrors in all directions as in the previous figures. 

From Fig. [5], we observe that at first the fixed bed floats like a single cluster, but the 
cluster becomes unstable. Then the fixed bed becomes a fluidized bed. In this fluidized 
bed, we observe the formation of bubbles periodically at the bottom, which float up, and 
travel through the bed. We also observe particle convection around bubbles. This kind of 
bubble formation from fixed beds is similar to that in the real experiments j57] and large 



size simulation fTl|]15] -|17| . 



Figure || shows the standard deviation of the particle velocities and Figure |7| shows the 
number of particles in the region at 8/25 of the cell height above the bottom. These figures 
show that bubble formation occurs periodically at the peaks of the standard deviation. 

At the end of this section, we summarize important characteristics in fluidized beds. 
When we compare the results in Figs. |2| and |] or [5], we observe that the particle motion 
has different characteristics. In Fig. |], particles float up in dilute regions and fall down in 
concentrated regions. This tendency can be understood as follows. The concentrated regions 
can be regarded as clusters, where flow cannot penetrate inside clusters. Thus clusters may 
have fast sedimentation velocity as in the definition of V oc a 2 with the radius a. On the 
other hand, slugging motion is understood from the standard theory of sedimentation of 
homogeneous suspensions, where dilute regions have larger sedimentation rate than dense 
regions. In Fig. ^ we see no isolated clusters, and particle distribution for horizontal 
direction is almost uniform. If particles are dispersed uniformly, fluid flow feels larger drags 
from many particles than that from small particles. In other word, the difference between the 
bubbling phase and the slugging phase is whether or not the convection of particles exists. 
To produce convections, we need both characteritics of clusterings and sedimentation. Thus, 
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the particle motions in real systems are determined by complex combinations of two different 
characteristics. 



IV. ANALYSIS 

In this section, we analyze the data obtained from our simulation. At first we discuss 
the velocity distribution functions (VDFs), where not only Gaussian-like but non-Gaussian 
distributions are observed and close relationship between the deviation from Gaussian and 
the local density of particles is indecated. Next we discuss the power spectra obtained from 
the data, which suggest the existence of tails obeying power laws. These power law tails 
may correspond to those reported by Taguchi |ft8|| , and taken as an evidence for powder 
turbulence |39| . 



We display VDFs in our simulation for systems with the square periodic cells in Fig. |8|. 
From this figure, it is obvious that the VDFs are near Gaussian but anisotropic. In this case, 
the vertical VDF has two branches, each obeying a different Gaussian distribution. This 
transmission in the VDF has also been observed in the simulation of the two-fluid model 
|19| . The anistropic poperty of VDF indicates that the system does not have any local 



equilibrium. 

For the systems with fixed particles, the form of VDFs is similar to an exponential one 
(Fig. §). It is interesting that the introduction of fixed particles produces an exponential 
VDF. In this system as shown in Fig. |6], active states which has large standard deviation 
of particle velocities and inactive states of motion of particles emerge in turn. From Figs. 
H and [7|, we indicate that the active state corresponds to the appearance of a bubble and 
in the inactive state contains is no bubble. To see the quantitative difference between two 
states, we investigate VDF in each state separately. Here we define the active state as the 
time region where the standard deviation of particle velocities is more than OAUq, and the 
inactive state as the period where the standard deviation is less than 0.3Uo- Figure |1(] is the 
VDFs for the active and the inactive state. From this figure, VDF in the inactive state is 
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close to an exponential distribution, while that in the active state is a Gaussian distribution. 
For the systems with slugs, VDFs are far from Gaussian futher and also different from 



simple exponential distribution (Fig. [11]). We attempt to fit them to t-distribution, 



f{U)~ll + aU*)~ b , (22) 



where a and b are parameters. This t-distribution is observed by Taguchi and Takayasu 
4Cj for VDFs in vibrating beds and by Sinai and Yakhot [f|lj for passive scalars in turbu- 



lence. From this figures, although the horizontail VDF can be fitted by the t-distribution, 
the vertical VDF is far from any t-distribution. However the tail for negative U y can be 
understood by the following simple picture. In the slugging state, most of particles are 
included in the dense region and only a few particles are falling in the slug. Here we assume 
that falling particles start from zero relative velocity to the dense region. The particles are 
accelerated by the gravity then they collide with the top of the dense region and join in it. If 
this assumptions are valid, the VDF of falling particles is uniform in the range between the 
initial and the terminal falling velocities, which may correspond to the tail for U y observed 
in Fig. [TTJ . To investigate the VDF in the dense region, we substract the uniform distri- 
bution from the original VDF (Fig. |TT] (b)) in the range of —0.4 < U < —0.115. Here the 
probability of the uniform distribution is approximated by the original VDF in the range of 
—0.7 < U < —0.4, and U = —0.115 is the location of the peak of the original VDF. The 
resultant VDF is shown in Fig. [12| in the range of —0.4 < U < 0.3. We can fit this well 
by (|22|) for each side of the peak, although there is still anisotopy in the figure. Thus, t— 
distribution seems to be applicable to our system. 

From our results of VDFs, we may indicate that non-Gaussian properties are closely 
related to the local density of particles. Because in relatively dilute case such as the system 
with square unit cells and the active state in the system with fixed particles, VDFs are 
close to a Gaussian, while in the dense case such as the inactive state in the system with 
fixed particles and the system which has slugs, VDFs are far from Gaussian. These two 
characteristics may be understood as follows. In active states, the density of particles are 
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relatively low. As a result, particles can move almost freely without influence of lubrication 
effects. Then collisions among particles occur at random, and inelastic effects from viscous 
terms which is mainly from the lubrication force may be suppressed. On the other hand, 
in inactive states, the density of particles is extremely high, and the lubrication effect also 
becomes important. Therefore, the inelastic effects dominate the randomness to produce 
Gaussian distributions. 

It is interesting that we observe non-Gaussian, exponential-like, VDFs in our systems. 
Non-Gaussian probability distributions are observed in various systems. In fluid turbulence 
non-Gaussian distributions are found in the probability distribution functions of velocity 



differences [[42] and passive scalars |43|,|44| . In simulations of vibrating beds, VDFs of particles 



are found to be described by a non-Gaussian VDFs | 40|| . For the simple models consisting 



of hard-spheres, VDFs also show non-Gaussian distribution |4qj4q] . In astronomy non- 
Gaussian VDFs are also reported . 



We also indicate that in the region obeying non-Gaussian VDFs, the kinetic temperature 
which is defined through the deviation of the Gaussian VDF cannot be used. Although it 
is possible to introduce the granular temperature from the different context, we need to be 
aware of diffrence between this granular temperature and the usual kinetic temperature. 

Next, we investigate the power spectra in frequency E{lo) defined through 

N 

N 



£ m4E<U»-U»>, (23) 



where 

U(u) = / dt e~ luJt V(t). (24) 



The results shown in Figs. O, [T3 and [T5| are obtained by a standard fast Fourier transform 



routine with the Parzen window fi8 |. All of these figures indicate that there are three 
regions; the spectra seems to be white in low u, there are some peaks in ther middle, and 
the spectrum obeys a power law in high u. To understand the mechanism to make three 
regions, we consider the following three characteristic frequencies, 
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uj l = 2nj, (25) 
u P = % (26) 

"r = ^ (27) 

where we adopt the average particle velocity relative to the fluid for the systems without 
fixed particles and the induced flow rate for the system with fixed particles as U. T p and 
T r introduced in ([?[) are the passing time of the particle scale and the relaxation time 
respectively. These frequencies are also shown in the figures. From these figures, ujl and u r 
seem to correspond to boundaries of three regions; the white spectrum region, the region 



with some peaks, and the region obeying power law. The peak near lo/luo ~ 13 in Fig. |TJ 
corresponds to the frequency of bubble formations, where uq = 27r/2048 x 10 _4 [sec]is the 
smallest frequency produced from our entire simulation. In the higher frequency range, all 
of these figures show the existence of a power law E(u) ~ uj 13 between u r and uj p whose 
exponent is fluctuated between —1.49 and —1.63. This is not far from the Kolmogorov 
spectrum E(u) ~ tu~ 5 / 3 in fluid turbulence. About the power laws, we will discuss futher in 
the next section. 

V. DISCUSSION 

Now we discuss the origin of the power laws observed in the power spectra, which are 



similar to Kolmogorov scaling [ 49| . Taguchi [^8j also observed the Kolmogorov-like scaling 
in vibrating beds. Therefore we need to clarify whether the power law observed in our 
simulations is the same as that of Taguchi. 

Our system contains an energy source on the scale of the lowest wave number, and the 
energy is dissipated in the fluid on the scale of the highest wave number. Vibrating beds have 
also common feature. Therefore, we might expect to have a cascade process, as assumed by 
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Kolmogorov and proposed by Richardson, where the energy dissipation rate e determines 
the statistical properties of small particle motion. For this purpose, we assume that the 
motion of particles is determined by only the energy dissipation rate e and the viscousity \i. 
As discussed in Sec. ||, the particle motion is described by only one relevant dimensionless 
parameter, the Stokes number, St = 2p p aV/9fi in the case of Pe ^> 1 and Re <C 1. Although 
the Stokes number may be regarded as the effective Reynolds number for the particle fluid, 
it contains only the energy dissipation length scale. Instead of St it is convenient to use the 
effective Reynolds number for particle fluid Re^ = VL/u p = 9LSt/2a defined in Eq.(8). 
From dimensional analysis, we obtain the following scaling for the energy spectrum 

E(k)=e^F(^-,Re {p) y (28) 

where F(x, y) is a dimensionless function, kd is the dissipation scale which is the maximum of 
a -1 and the Kolmogorov scale e 1//4 z/~ 3 / 4 . In the region of high wave numbers we might guess 
that there is a balance between energy injection and cascade process, where the dissipation 
is not important. Therefore, we obtain 

E(k) = e 2/3 k~ 5/3 . (29) 



This result is equivalent to that by Kolmogorov [[49[ . 



In order to ensure the scaling ansatz, we need to impose the condition 

Lk d ~ Re 3 ^ > 1. (30) 

Equation ([30]) means that the inertia spectrum can be observed in our system (Table ||). 
Although there is an ambiguity in choosing L, it is interesting that the values of -Re^ 4 
correspond to the strength of the non-Gaussain VDFs. We also expect the cutoff a -1 is 
larger than the dissipation scale e 1//4 ^~ 3//4 , because we do not observed any dissipation range 
(see also [[38]). We also comment on the reason to observe turbulent charactersitics in 
relatively small Re^, where in usual fluid, the value of Re^ in Table |TJ is not enough 
large to show turbulent charcateristics. We guess that inelatic scatterings of particles in our 
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system enhance chaotic characteristics, while the molecules in pure fluid only have elastic 
scatterings. We can use this argument to vibrating beds, where /i is replaced by the effective 
viscosity of particle flow. 

We only discuss the spectrum for E(k) and not for our observed E(u). Although the 



for the inertia range in both Euclerian and Lagrangian coordinate frames can be derived 
from (p9|) with the assuption of space-time symmetry. His argument is independent of the 
basic equations and assumes that the system can be charcaterized only through the kinetic 
viscoisty and e as in our argument. Therefore, we may expect that his argument is applicable 
to our case. 

We also comment on the reason that we observe the Kolmogorov like scaling even in 
two-dimensional systems. In our case the enstrophy is not a conserved quantity in two 
dimensional case. Therefore it is not surprizing that we observe the Kolmogorov-like scaling 
in two-dimensional systems, because the derivation of does not contain any information 
about the spatial dimensionality. 



In this paper, we have succeeded in simulating granular systems with fluid flow without 
any empirical parameters. Our simulation reproduces realistic slugs and bubbles. According 
to the results of our simulation, we confirm that particle convection is important for bubble 
formation. From our simulation we observe several velocity distribution functions, Gaussian 
and non-Gaussian or exponential distribution. We also observe power spectra obeying power 
laws. We briefly explain the relationship between our observed spectra and Kolmogorov 
scaling. In this sense, our paper has confirmed the universality of the powder turbulence 
proposed by Taguchi P5fl . 



general relationship between E(k) and E(u) is unclear, Viecelli [|5tJ demonstrates 




(31) 



VI. CONCLUSION 
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Before closing this paper, we comment on our approximation in using the Stokes equation 
for fluid flow. In real experiments and simulations by engineers, the particle radius, typically 
a ~ 10 _1 [cm], is often much larger than that we assume (a ~ 10~ 3 [cm]). This is because 
the direct interaction among particles except for hard-core interactions dominates all other 
forces, and it is difficult to observe fluidized beds in such small particles For the 

fluidized beds with such large particles, we need to consider the advection term in the 
fluid motion and the drag which should contain a term in proportion to the square of the 
difference of velocities between fluid and particles ||13|| . To simulate this system without 
introduction of empirical laws is almost impossible. We, however, have shown that these 
kinds of complexities are irrelevant for the formation of bubbles and slugs. The essence of 
the physics in fluidization can be understood when we drop irrelevant complexities such as 
direct interparticle forces and the advection term in fluid flows. 
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APPENDIX A: HYDRODYNAMIC INTERACTION 

In this section, we explain how to treat hydrodynamic interactions. When we adopt 
the Stokes approximation, the calculation of the hydrodynamic interactions is equivalent to 
constructing the resistance matrix. Our aim is to calculate the resistance matrix with ac- 
ceptable accuracy taking into account computational efficiency We thus adopt the Stokesian 
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dynamics developed by Brady and his coworkers |26|^28| to describe colloidal dispersions. 
They distinguish hydrodynamic long-range interactions from the lubrication force which 
represents short range hydrodynamic repulsive interactions. For the long-range part, we 
use the multipole expansions, while for the lubrication part, we use the exact solution of 
two-body problem pTI, 



as 







XJW - u°° 




— R-2B " 




p( 2 ) 




U( 2 ) - U°° 



(Al) 



where u°° is the fluid velocity without particles. The lubrication part is separated from the 
exact solution as 



R l 2B - R 2B - (M^) 



-1 



(A2) 



where is the Rotne-Prager tensor ]5B[ , which represents the long-range interaction. Let 
R lub be the linear combination of all possible pairs of R^. Then, an approximate resistance 
matrix is represented by 

R = (M 00 )- 1 + R lub . (A3) 
The validity of this method has been confirmed in various numerical experiments for colloid 



systems. This idea is also valid in the theoretical calculation of the sedimentation rate p5 

For the long-range interaction, we apply periodic boundary conditions to describe the 
system containing infinite number of particles. For simplicity, we consider only the con- 
tribution of force. Note that higher order corrections from torque have been discussed by 



Brady et al. [28]. Beenakker |54[] has shown that the mobility matrix can be written Ewald's 
summation as 



N 



6iifia(Ui 



U; 



7 /9=1 
N 



1 

y EE «* (k A ■ (x« - ^)) Mg\k x )Ff - Mg\v = 0)F/ 



(A4) 



A^0/3=l 

where N is the number of particles in the unit cell and M^\r), M^^k) and M^\r = 0) 
are respectively given by 
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, ,<T)/ x r , f r /3a la 3 \ _ /3 a 3 a 3 " 
M«(r) = erfc(er) U 3 (_ + -_] + (-- - -- 



7T 



(4a 3 <fr 4 + 3a£ 3 r 2 - 20a 3 £ 5 r 2 - ^a£ + 14a 3 £ 3 + « 3 ^) 



/ 3 i \ 

+f i f i ( -4a 3 «fr 4 - 3a£ 3 r 2 + 16a 3 £ 5 r 2 -a£ - 2a 3 £ 3 - 3a 3 £— J 



(A5) 



M§\k) = 6vra (<f y - kh) y 2 (l - f 1 + ^ + ^) e ~^' ( A6 ) 



and 



Af<?> (r = 0) = ^| ( 6 a^ - f a 3 e 3 ) . (AT) 



Here erfc(x) is the complimentary error function given by 

2 r 00 



erfc(x) = —j= / exp(— z 2 )dz, (A8) 

\/7T Jx 



and £ with units of inverse length is an arbitrary parameter which we choose to minimize the 
number of lattice-sums in both real space and k-space. In the simulation, we use £ = 
where L is the average of the length of the periodic cell. The suffix 7 = (ni,n 2 , n 3 ) represents 
a periodic cell in the real space and r 7 is the lattice vector given by 

r 7 = (n 1 Li,n 2 L 2 , n 3 L 3 ) , (A9) 

where L\,L 2 , L 3 denotes the length of unit cell in each direction. The wave vector k\ in the 
reciprocal cell A = (m 1 ,m 2 ,m 3 ) is given by 

27rmi 2 r nm 2 2iTm 3 



(A10) 



The summation M-j in Eq. (|A4|) means that the contribution of the sum in which a = (3 
at 7 = is eliminated. The contribution of k = in the lattice sum in reciprocal space is 



canceled by that of the average force acting on the fluid [^3 • 

We have checked the validity of our program to examine several sedimentation velocities 
under regular configurations like simple cubic, body-center cubic and face-center cubic. We 
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compared our results with the Stokesian dynamics by Brady et al. p8| in which we choose the 
corresponding model, neglect higher order moment of hydrodynamic interactions, and with 
the exact solution by Zick and Homsy |55[|. From this test, we recover the corresponding 



result by Brady et al. ||28|| , that is, our program seems to work correctly. Furthermore, 
our result is very close to the exact solution by Zick and Homsy ]55] except for extremely 
concentrated region. Thus the simplification of neglecting higher order moments does not 
cause serious problems. 



APPENDIX B: FIXED PARTICLES 

In this section, we show how to calculate the terminal velocities when we introduce 
particles fixed in the space. If we obtain the terminal velocities, the motions of free particles 
are determined by fl2~I|) . 

If we know the resistance matrix which contains both free particles and fixed particles, 
we obtain the following equation, 



F 








u m 


- u°° 

m 


_F/_ 




R/m R// 






— 11 00 



(Bl) 



where the suffix m represents free particles and the suffix / represents fixed particles. We 
already know the force acting on free particles F m , the velocity of fixed particles U/ = 
and the velocity of induced fluid u°°. ( |B~1| ) can be solved for U m as 

U m - iC = R~^ ■ [F m + R mf ■ uf^j . (B2) 

Therefore the terminal velocities of free particles can ben represented by known variables. 
This procedure is applicable to systems with periodic boundary conditions. In fact, from 
Eq. flAH ) the mobility matrix can be composed of particles in the unit cell. 
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FIGURES 

FIG. 1. The relationship between particle radius and dimensionless parameters. The straight 
lines connects the points in which dimensionless parameters are 1. 

FIG. 2. Snapshots of the slug. The number of particles in the unit cell is 72 and the volume 
fraction is 0.348. The ratio of the height to the width is 6. Time interval is 2.0 x 10~ 2 sec. 

FIG. 3. Snapshots of the slug in three-dimensional simulation. The number of particles is 50 
in the unit cell and the volume fraction is 0.45. The ratio of the height to the width is 4, and time 
interval between them is 2.0 x 10~ 2 [sec]. 

FIG. 4. A snapshot of the simulation with square cells with the periodic boundary condition. 
It is shown with four periodic images. The number of particles in the unit cell is 90 and the volume 
fraction is 0.327. 

FIG. 5. The sequence of the simulation with fixed particles. The number of free particles 
and fixed particles in the unit cell are 128 and 10, respectively. The velocity of induced fluid is 
u°° = 0.3Uo, where Uq = mg/67Tfj,a is the one particle sedimentation velocity. The time interval 
from the left-up to the right-down is 0.112 sec. We can see three bubbles flowing up through the 
bed (except for the initial instabilities). 

FIG. 6. The time evolution of standard deviation of particle velocities in the simulation of 
Fig. |. Velocit ies are normalized by one particle sedimentation velocity Uq. 1 step means 10 sec. 

FIG. 7. The time evolution of the volume fraction of particles in the area that is the horizontal 
region at 8/25 of the cell height above the bottom in the simulation of Fig. ||. 1 step means 10~ 4 
seconds. 
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FIG. 8. The velocity distribution functions (VDFs) of (a) the horizontal direction and (b) 
the vertical direction in the simulation of Fig. || (with square unit cells). The velocities are 
normalized by one particle sedimentation velocity Uq. The vertical component of velocity is plotted 
in the coordinate system in which the average velocity of the fluid is zero and the direction of 
gravity is negative direction. The horizontal VDF can be approximated by exp(— 20.0U 2 ), where 
U = U /Uq. The left hand side and the right hand side of the vertical VDF can be approximated 
by exp(-10.0(*7 + 0.2) 2 ) and exp(-4.5(J7 + 0.2) 2 ), respectively. The VDF f(U) is normalized by 
JdU f(U) = l. 

FIG. 9. The velocity distribution functions (VDFs) in the simulation of Fig. ||(with fixed 
particles). The velocities are normalized by one particle sedimentation velocity Uq and we use the 
coordinate system in which the fixed particle's velocity become zero. The VDFs in horizontal di- 
rection and vertical direction can be approximated by exp(— 5.0|?7|) and exp(— 3.6|J7|), respectively. 
Here U = U /Uq. The normalization is the same as that in Fig. ||. 

FIG. 10. The velocity distribution function (VDF) in (a) the active state and (b) the inactive 
state of the simulation with fixed particles which is the same of that in Fig. [9|. The velocities are 
normalized by one particle sedimentation velocity Uq and we use the coordinate system in which 
the fixed particle's velocity become zero. The Gaussian-fitted function is also ploted. 

FIG. 11. The velocity distribution function (VDF) of (a) the horizontal direction and (b) 
the vertical direction in the simulation of Fig. ^ (with rectangular unit cells). The velocities 
are normalized by one particle sedimentation velocity Uq. The VDF in the horizontal direction 
can be approximated by exp(— 25.0|C/|) in the exponential form or (1.0 + 700. 0?7 2 ) -1 ' 5 in the 
t-distribution form. The VDF in the vertical direction for the right-hand-side by can be approx- 
imated by exp(-50.0|J7 + 0.115|), in the exponential form or by (1.0 + 800.0(£7 + 0.115) 2 ) -20 in 
the t-distribution form. Here U = U /Uq. The normalization is the same as that in Fig. ||. 
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FIG. 12. The substracted velocity distribution function (VDF) from the VDF in Fig. 11 
(b). We make this figure by substracting the uniform distribution, representing the contribution of 
falling particles in the slug, from the original VDF in the range of —0.4 < U < —0.115. This VDF 
can be fitted for the left-hand-side by (1.0 + 700.0(^7 + 0.115) 2 ) -1 ' 5 and for the right-hand-side by 
(1.0 + 800.0(L> + 0.115) 2 ) -2 ' . Here U = U/U . 

FIG. 13. Power spectrum E(u) of (a) the holizontal velocity and (b) the vertical ve- 
locity in the simulation of Fig. |4| (with square unit cells). The frequency is normalized by 
ujq = 2-7r/2048 [step -1 ] ~ 30.7[sec -1 ]. The straight lines are least-square fits between the frequency 
from u> r to top with —1.597 ± 0.009 and —1.626 ± 0.009 slope for U x and U y respectively. 

FIG. 14. Power spectrum E{uj) of (a) the holizontal velocity and (b) the vertical veloc- 
ity in the simulation of Fig. ^ (with rectangular unit cells). The frequency is normalized by 
ujq = 2-7r/2048 [step -1 ] ~ 30.7[sec -1 ]. The straight lines are least-square fits between the frequency 
from u> r to top with —1.494 ± 0.01 and —1.513 ± 0.009 slope for U x and U y respectively. 

FIG. 15. Power spectrum E(uj) of (a) the holizontal velocity and (b) the vertical ve- 
locity in the simulation of Fig. |5|(with fixed particles). The frequency is normalized by 
loq = 2-7r/2048 [step -1 ] ~ 30.7[sec -1 ]. The straight lines are least-square fits between the frequency 
from uj r to top with —1.514 ± 0.007 and —1.490 ± 0.007 slope for U x and U y respectively. 
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TABLES 

TABLE I. Values of the dimensionless parameters and the characteristic time T r [sec] and 
T p [sec] when the particle radius is 10 _3 [cm]. 
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TABLE II. Values of Re^ in our simulations, where L is measured by the maximum length 
of the unit cell. 
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